Dynamics of simulated water under pressure 
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We present molecular dynamics simulations of the SPC/E model of water to probe the dynamic 
properties at temperatures from 350 K down to 190 K and pressures from 2.5 GPa (25 kbar) down 
to -300 MPa (—3 kbar). We compare our results with those obtained experimentally, both of which 
show a diffusivity maximum as a function of pressure. We find that our simulation results are 
consistent with the predictions of the mode-coupling theory (MOT) for the dynamics of weakly 
supercooled liquids - strongly supporting the hypothesis that the apparent divergences of dynamic 
properties observed experimentally may be independent of a possible thermodynamic singularity at 
low temperature. The dramatic change in water's dynamic and structural properties as a function 
of pressure allows us to confirm the predictions of MCT over a much broader range of the von 
Schweidler exponent values than has been studied for simple atomic liquids. We also show how 
structural changes are refiected in the wave-vector dependence of dynamic properties of the liquid 
along a path of nearly constant diffusivity. For temperatures below the crossover temperature of 
MCT (where the predictions of MCT are expected to fail) , we find tentative evidence for a crossover 
of the temperature dependence of the diffusivity from power-law to Arrhenius behavior, with an 
activation energy typical of a strong liquid. 



I. INTRODUCTION 

The "slow dynamics" and glass transition of both sim- 
ple and molecular liquids has been a topic of significant 
interest in recent years. The initial slowing down of liq- 
uids at temperatures down to Tc « \.2Tg, where relax- 
ation times approach 1 ns, has been well described by the 
mode-coupling theory (MCT) |l|. MCT has been suc- 
cessfully applied to a wide variety of real and model sys- 
tems including hard spheres Ni8oP2o @|, Si02 jg, 
and polymer melts However, there has not been an 
extensive test of the validity of the MCT predictions for 
a model system over a wide range of pressures and along 
different thermodynamic paths. 

At low pressure, it was shown previously that the 
power-law behavior of dynamic properties in the SPC/E 
model can be explained using MCT Furthermore, 
the possible relationship between the experimentally- 
observed power-law behavior and the predictions of MCT 
has been discussed ||^-[ll|]. The experimentally-observed 
locus of apparent power-law singularities of dynamic and 
thermodynamic properties [Fig. |l|] is of particular inter- 
est [|l^,|l^, and has catalysed the development of three 
scenarios to explain the anomalous properties of water: 
(i) the existence of a spinodal bounding the stability of 
the liquid in the superheated, stretched, and supercooled 
states p2| , p^ ; (ii) the existence of a liquid- liquid phase 
transition line separating two liquid phases differing in 
density [^-|l9|; (iii) a singularity- free scenario in which 
the thermodynamic anomalies are related to the pres- 



ence of low-density and low-entropy structural hetero- 
geneities 1 20 . The predictions of MCT are of interest 



since MCT might account for the apparent power-law 
behavior of dynamic properties on cooling, thereby re- 
moving the need for a thermodynamic explanation of the 
dynamic properties of water. 



In this article we focus on two related issues: (i) the 
possibility of using MCT to explain the slow dynamics 
of water under pressure, and (ii) a test of the validity 
of the MCT predictions over an extremely wide pressure 
range in a system with dramatic structural changes. We 
find that MCT provides a good account of the slow dy- 
namics of the SPC/E model for water at all pressures, 
with the structure evolving continuously from an open 
tetrahedral network to a densely packed fluid, similar to 
a Lennard- Jones type liquid. By examining the wave- 
vector dependence of collective dynamics, we are able 
to discover how these structural changes are reflected in 
the dynamic behavior of the liquid. We are also able to 
test the validity of the relationship predicted by MCT for 
the diffusivity exponent 7 and the von Schweidler expo- 
nent 6 over a wide range of values 7 and b [Fig. j^]. Our 
results support the predicted relationship of these expo- 
nents. A brief report of a subset of the present results 
for the SPC/E potential has recently appeared [^l[|. The 
dynamic properties of the ST2 potential [^ at one pres- 
sure, and of the TIP4P potential at several pressures, 
have also recently been discussed. 
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FIG. 1. Phase diagram of water. The extrapolated diver- 
gence of the isothermal compressibility (o) and the ex- 
trapolated divergence of D {•) The different loci of these 
two singularity lines are consistent with the possibility that 
the two phenomena may arise from different explanations. 
Also shown are the melting line (Tm) and coexistence lines of 
several ice polymorphs and the experimental limit of super- 
cooling (Th). 



II. MODE COUPLING THEORY 

We will focus our discussion on the idealized form 
of MCT, originally formulated to describe spherically- 
symmetric potentials. Recent extensions have been made 
to account for the rotational motion present in non- 
spherical molecular systems [|4|, such as water. The ide- 
alized version of MCT has been shown to provide a good 
account for the center-of-mass motion for the SPC/E 
model l^,^. We provide only a brief account of the 
MCT predictions relevant to the results of this article, 
and we refer the reader to extensive reviews for more 
information §,|j2^. 

MCT assumes that localization, or "caging" of 
molecules due to the slow rearrangement of neighboring 
molecules, is the source of the dramatic increase of re- 
laxation times on cooling, leading to a strong coupling 
between single particle motion and the density fluctua- 
tions of the liquid. Indeed, according to MCT, the static 
density fluctuations, measured by the structure factor 
S{q)^ entirely determine the long time dynamic behavior. 
MCT accounts for the loss of correlation by the interac- 
tion of density mode fluctuations, ignoring other possible 
mechanisms for relaxation. MCT predicts the asymptotic 
power-law divergence of correlation times, and power-law 
vanishing of the diffusion constant 

D Do{T/T,-l)'' (1) 

at a critical temperature Tc — Tc{P), where we refer to 
7 — 7(P) as the diffusivity exponent. In real systems, 
"freezing" of the system dynamics is avoided at Tc, as 
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FIG. 2. The line shows the predicted relationship between 
b and 7 from MCT. The symbols show the calculated values 
for the SPC/E model: (o) from this work, (filled O) from 
Ref. §. 

relaxation mechanisms not accounted for by MCT be- 
come significant. However, Tc can still be interpreted as 
a "crossover temperature" where the dynamics change 
from being dominated by density fluctuations to being 
controlled by "activated" processes. Some recent work 
has also demonstrated the significance of Tc as a crossover 
temperature where relaxation occurs primarily through 
basin hopping [P7|-p0|, in the energy landscape view of 
liquid dynamics IP 5[p2| . 

MCT predicts that the Fourier transform of the 
density-density correlation function or intermediate 
scattering function 

F(g,O^^^E e-'^-I-W-^Wl^ (2) 

decays via a two-step process. In the first relaxation step, 
F{q,t) approaches a plateau value FpiatoauC?) which is 
described, to leading order in time, by a power law with 
exponent a — a(P), 

F(g,t)-Fpiatcau(9)-^i", (3) 

At larger times, F{q,t) decreases from Fpiatcau(<7) and 
MCT predicts the decay obeys the von Schweidler power 
law to leading order in time 

Pplatcau(g) -P(<?,i) (4) 

where b = b{P) is known as the von Schweidler exponent. 
The region of validity of Eqs. (|) and (|) can be quite 
limited. 

The slow relaxation of F{q, t) has a characteristic re- 
laxation time r that is also predicted to have asymptotic 
power law dependence on temperature. 
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r~To(T/T,-l) 



(5) 



with the same vahie of the exponent 7 as for the diffu- 
sion constant. Hence, Eqs. (|l|) and (|^) predict that the 
product Dt is not singular as T — > Tc, hence we take the 
product to be constant over the range that Eqs. (|l|) and 
are vahd (neglecting corrections to scaling). 
MCT predicts that the scaling exponents a, h, and 7 
are not independent; a and b are related by the exponent 
parameter A using the relationship 



[r(i-a)p [r(i-6)]= 



r(l-2a) r(l + 26) 



(6) 



where T{x) is the gamma function. MCT also relates 7 
to a and 6 via 
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Because of Eqs. (||) and (|^), only one exponent value 
is needed to determine all others, so calculation of two 
exponents determines if the dynamics of a system are 
consistent with the predictions of MCT. Furthermore, 
these exponents are expected to depend on the path along 
which Tc is approached. 

After F{q,t) departs from the plateau, F{q,t) is well- 
described by a Kohlrausch- Williams- Watts stretched ex- 
ponential 



F{q, t) ^ A{q) exp 



t 



riq) 



(8) 



where r(q) is the relevant relaxation time. Moreover, it 
has been shown that the exponent (3 — (3{q) is related to 
the von Schweidler exponent M 



lim (3{q) = b. 



(9) 



This relation facilitates evaluation of b, since the region 
of validity of Eq. fl) is difficult to identify in practice. 



III. SIMULATIONS 

We perform MD simulations of 216 water molecules in- 
teracting via the SPC/E pair potential 0. The SPC/E 
model treats water as a rigid molecule consisting of three 
point charges located at the atomic centers of the oxy- 
gen and hydrogen, which have an OH distance of 1.0 A 
and HOH angle of 109.47°, the tetrahedral angle. Each 
hydrogen has charge qn — 0.4238e, where e is the fun- 
damental unit of charge, and the oxygen has charge 
qo = ~2qH- In addition, the oxygen atoms of separate 
molecules interact via a Lennard- Jones potential with pa- 
rameters a = 3.166 A and e = 0.6502 kJ/mol. 

Our simulation results are summarized in Table |. For 
T < 300 K, we simulate two independent systems to im- 
prove statistics, as the long relaxation time makes time 



averaging more difficult. We equilibrate all simulated 
state points to constant T and p by monitoring the pres- 
sure and internal energy. We control the temperature us- 
ing the Berendsen method of rescaling the velocities , 
while the reaction field technique with a cutoff of 0.79 
nni accounts for the long-range Coulonibic interac- 
tions. The equations of motion evolve using the SHAKE 
algorithm |37) with a time step of 1 fs, except at T = 190, 
where a time step of 2 fs is used due to the extremely slow 
motion of the molecules. Equilibration times at high tem- 
peratures are relatively small. At low T, extremely long 
equilibration times are needed. The structural and ther- 
modynamics properties may be obtained after relatively 
short equilibration times. However, dynamic properties 
show significant aging effects (i.e. dependence of mea- 
sured properties on the chosen starting time) if great care 
is not taken in equilibration. 

For production runs, it is desirable to make mea- 
surements in the isoenergetic/isochoric ensemble (NVE). 
However, a small energy drift is unavoidable for the long 
runs presented here, so we again employ the heat bath of 
Berendsen, using a relaxation time of 200ps |^ . The large 
relaxation time prevents an energy drift but achieves re- 
sults that are very close to those that would be found 
if it were possible to perform a simulation in the NVE 
ensemble. 

Since we perform long runs for many state points, we 
store the molecular trajectories {ri,pi} at logarithmic 
intervals to avoid storage problems that linear sampling 
presents. Specifically, we sample configurations at times 
growing in powers of 2 up a maximum time tmax- We be- 
gin a new sampling cycle each time imax (relative to the 
cycle starting time) is reached. This sampling method al- 
lows for calculation of dynamic properties on time scales 
spanning eight orders of magnitude (from 1 fs to 100 ns) 
using a relatively small amount of disk space. Still, more 
then 2 GB of storage was required for storing configu- 
rations at T = 210 K. Our simulations have a speed 
of approximately 200 /zs per update per molecule on a 
MIPS RIOOOO processor, representing a total calculation 
of approximately 8.4 years of CPU time, including the 
systems of 1728 molecules discussed in the appendix. For 
the larger systems, we utilize a parallelized version of our 
simulation code on eight processors to improve perfor- 
mance. 



IV. STATIC STRUCTURE FACTOR 



We first summarize the structural properties of our 
simulations in order to better understand the relation- 
ship between the changes in structure with the changes 



in dynamic behavior, which we will detail in Sec. VIII 



Other studies have considered the structural and ther- 
modynamic properties of SPC/E in a large region of the 
(P.T) plane ]38|-^, so the present discussion is brief. 
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FIG. 3. The oxygen-oxygen structure factor S{q): (a) De- 
pendence on p for T = 210 K. (b) T-dependence along the 
p = 1.0 g/cm^ isochore. Notice that changing T has little 
effect, while changing p has a more pronounced effect. 



The MCT theory requires as input the static density- 
density correlation functions. In the case of water, the 
structure of the system is very sensitive to the value of 
the external control parameter {P, T) . Hence, for all state 
points simulated, we calculate the oxygen-oxygen partial 
structure factor lISSll 



S{q) 



1 



N 



(10) 



Several studies have carefully calculated the structure of 
simulated water, and found surprisingly good agreement 
with experiments At T = 210 K, we show the 

structural changes from low to high density [Fig. ||(a)]. 
The structure at low density /pressure is similar to that 



observed for low-density amorphous (LDA) solid water, 
consisting of an open tetrahedral network. At high den- 
sity/pressure, water is very similar to high density amor- 
phous (HDA) solid water, where core-repulsion domi- 
nates, similar to simple liquids under pressure. 

We show the evolution of S{q) as a function of T along 
the p = 1.0 g/cm'^ isochore in Fig. ^(b). We note that 
in the temperature range from 190 K to 300 K, where 
the dynamics show the most dramatic change in behav- 
ior, S{q) shows only small changes in the first two peaks. 
Also, the location of the first maximum qq in S{q), the 
wave vector at which F(q, t) typically shows the slowest 
relaxation, does not appear to change significantly. All 
other densities and temperatures show a relative smooth 
interpolation of Figs. 0(a) and (b). 



V. MEAN-SQUARED DISPLACEMENT AND 
DIFFUSION 

The mean-squared displacement {r'^{t)) = (|r(t) — 
r(G)p) is shown in Fig. ^. All the curves show de- 
pendence at small time, as expected in the "ballistic" 
regime. For low T (e.g. T = 210 K [Fig. |(a)]), {r^{t)) 
shows relatively flat behavior over 3-4 decades in time. 
This is the "cage" region, in which a molecule is trapped 
by its neighbors and cannot diffuse, and is only vibrat- 
ing within its "cage." At low P, the cage consists of 
hydrogen-bonded neighbors in a tetrahedral configura- 
tion. This cage is relatively strong, compared to simple 
liquids, because of the H bonds. The size of the cage 
may be estimated by the value {r'^{t)) at the plateau, 
as shown in the inset of [Fig. ^a)]. Surprisingly, the 
size of the cage is not monotonic with density, and has a 
maximum at p « 1.1 g/cm^. We shall see that this corre- 
sponds roughly to the p at which D also has a maximum. 
We observe a small bump in {r'^{t)) at i « 0.35 ps, as ob- 
served in Ref . [^ . A system size study indicates that this 
may be attributed to finite size effects 

For long times, all the curves show linear t depen- 
dence, indicating that our simulations are in the diffusive 
regime. We extract the diffusion constant D using the 
asymptotic relation {r^(t)) = 6Dt. We plot the density 
dependence of D in Fig. |^ and find that the SPC/E po- 
tential, like water, shows an anomalous increase in D on 
increasing density. We also point out the feature that D 
shows a slight increase at very low density; namely, at 
p = 0.90 g/cm^ and T = 210 K. This can be attributed 
to the fact that the liquid is extremely stretched at this 
density, causing an increase in the defects of the bond 
network, and thus increased diffusivity. We show D as 
a function of pressure along several isotherms to com- 
pare with experimental measurements [Fig. ^ [^. The 
anomalous increase in D is qualitatively reproduced by 
our calculations for the SPC/E model, but the quan- 
titative increase of D is significantly larger than that 
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FIG. 4. Mean-squared displacement {r'^it)) for (a) all den- 
sities at r = 210 K and (b) all T along the p = 1.0 g/cm^ 
isochore. The inset of (a) shows the density dependence of 
the cage size. 

observed experimentally. This discrepancy may arise 
from the fact that the SPC/E potential is under- 
structured relative to water , so applying pressure al- 
lows for more bond breaking and thus greater diffusivity 
than observed experimentally. We also find that the pres- 
sure where D begins to decrease with pressure — normal 
behavior for a liquid — is larger than that observed ex- 
perimentally . This comparison of D with experiment 
leads us to expect that while the qualitative dynamic fea- 
tures we observe in the SPC/E potential may aid in the 
understanding of the dynamics of water under pressure, 
they will likely not be quantitatively accurate. 

We estimate D along the isobars P = —80 MPa, 
MPa, 100 MPa, 200 MPa, 300 MPa, and 400 MPa from 
the isochoric data. We confirm that along the -80 MPa 
isobar, our estimates agree with the —80 MPa calcu- 
lations of Ref. which employs the same truncation 
of the potential used here (see Sec. III). Along the MPa 
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FIG. 5. (a) Diffusion constant D along isotherms for each 
density simulated, (b) Relaxation time r of F{qo,t) along 
isotherms for each density simulated, (c) Test of the MCT 
prediction that Dt is constant along isochores. 
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FIG. 6. D as a function of pressure for various tempera- 
tures from (a) our simulations and (b) NMR studies of wa- 
ter i. 
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FIG. 7. Fit of each isochore to the power 
D ~ {TjTa - 1)^ predicted by MCT. We include the 
of Ref. |43] along the p = 1.0 g/cm^ isochore. 
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isobar, our estimates of D are smaller than those calcu- 
lated for SPC/E in Ref. |§, perhaps because Ref. |3|] 
chooses a different truncation of the electrostatic terms — 
highlighting the extreme sensitivity of the dynamics to 
changes in the potential. 

We fit D by the power law of Eq. (|]) along both iso- 
chores and the estimated isobars for T < 300 K [Figs. |^ 
and p. The values of the two fit parameters Tc and 
7 are given m Table @. We also include the data 
from Ref. Jp3[ along the p = 1.0 g/cm'^ isochore and 
from Ref. j^along the P = — 80 MPa isobar to improve 
the quality of the fits. At p = 1.40 g/cm^, we exclude 
T = 210 K when fitting D and obtain Tc = 209.3, since 
we expect the power law of Eq. to fail for T < Tc+5 K, 
because activated processes - such as "hopping" not ac- 
counted for in the idealized MCT - become significant 
and aid diffusion. To demonstrate the presence of hop- 
ping at p = 1.40 and T = 210 K, we plot the "self" 
part of the van Hove correlation function Gs{r,t), which 
measures the distribution of particle displacements r at 
time <, for several densities at T = 210 K [Fig. ^. For 
these densities where a power law adequately describes 
D, there is a single peak. At p = 1.40 g/cm^, we see 
a "shoulder" in Gs{r,t) at r « 0.2 nm in addition to a 
well-defined peak at r w 0.05 nm, indicating that particle 
hopping is significant. 

Along the p — 1.00 g/cm^ we have simulated to sig- 
nificantly lower T, allowing us to study the temperature 
dependence of D for T < Tc- Fig. |l^ shows that the low- 
est temperatures are consistent with the Arrhenius form 



(11) 




Arrhenius temperature dependence of D is not sur- 
prising in this region, since for T < Tc it is expected 
that the energy barriers the system must overcome to 



FIG. 8. Fit of the diffusion constant along each isobar to 
the power law D ~ (T/Tc-l)^ predicted by MCT. We include 
the data Ref. M along the -80 MPa isobar. 



rearrange exceed the thermal energy p7[ . Hence the mo- 
tion of the system is dominated by activated jumps over 
the energy barriers, as described by Goldstein [M. We 
obtain an activation energy of _E « 65 kJ/mol and ex- 
trapolate a glass transition temperature Tg « 125 K [Q, 
surprisingly close to the experimental value of 136 K. The 
extrapolated value of Tg is similar to that estimated in 
Ref. [4^ which studied hydrogen bond dynamics. More- 
over, our results are consistent with a crossover from 
"fragile" behavior (the behavior described by MCT) for 
T ^ Tc, to "strong" behavior (Arrhenius behavior with 
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FIG. 9. The van Hove correlation function G's(r, t) for sev- 
eral densities at T = 210 K. For each curve, t is chosen 
such that {r^{t)) ~ 0.1 nm^, well inside the diffusive regime 
(i.e., where {r^{t)) is linear in t). The presence of a pro- 
nounced shoulder in Gs{r,t) for p = 1.4 g/cm? indicates that 
hopping phenomena are significant, and thus deviations from 
power-law dependence are expected. 




FIG. 10. Arrhenius plot of D shows that the power-law 
behavior (solid line) appears to make a smooth crossover to 
Arrhenius behavior (dashed line) for T < Tc with activation 
energy _E « 65 kj/mol. This behavior may be related to a 
possible fragile-to-strong transition of the dynamic properties 
(see discussion in text). 

E w kBTg/25 « 40 kJ/mol for our estimate of Tg) for 
T ^Tc- The possibility of a "fragile-to-strong" crossover 
in water has been discussed recently based on experimen- 
tal findings ]4^ , but lower temperatures are required to 
test this possibility in the SPC/E model. 



VI. ISOCHRONES OF D AND THE LOCUS OF 

TciP) 



To construct isochrones of D (lines of constant D), 
we first estimate T{D) at values of D = 10~^cm^/s, 
2^0-5.5 cm^/s, lO^^cm^/s, and 10^''cm^/s, using the fits 
of Figs. 1^ and ^. Along the isobaric paths, we know al- 
ready P for these points, and along isochores we may 
estimate the value of P using the results presented in 
Table ||. We plot the isochrones in Fig. 

We also show the the loci of Tc{P) in Fig. |ll|(a), ob- 
tained from the fits in the previous section. We know P 
at Tc along the isobaric paths, and we estimate the P at 
Tc along isochores by extrapolating P in Table || to Tc- 
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FIG. 11. (a) Isochrones of D from simulation. The 
lines may be identified at follows: D = 10~^cm^/s (o); 
D = 10"^-^ cmVs (□); D = lO'^'crnVs (o); D = lO^^cmVs 
(A). The diffusion is also fit to D ~ (T - . The locus of 
Tc is indicated by (x). (b) Isochrones of D constructed from 
the experimental data in Ref. M. 
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FIG. 12. Pressure dependence of the difFusivity exponent 
7 defined by 7? ~ (T — TcY' ■ The symbols are as follows: 
(o) 7 measured from simulation along isochores; (□) 7 mea- 
sured from simulation along isobaric paths, which are esti- 
mated from the isochoric data; (filled O) 7 measured along 
the -80 MPa isobar in Ref. (filled A) experimental mea- 
surements of 7 in water from Ref. It is clear that the 
SPC/E potential fails to reproduce the qualitative behavior 
of 7 under pressure in liquid water. 



Using the experimental diffusion data of Ref. 
we also construct the behavior of the experimental 
isochrones following the same technique [Fig. |n|(b)] . The 
shape of the locus of TAP) compares well with that ob- 
served experimentally P, and changes slope at roughly 
the same pressure [Fig. ]ll| . Therefore, an explanation of 
the SPC/E dynamics using the MCT would support us- 
ing the MCT framework as an interpretation of the exper- 
imentally found locus of Tc{P). We find, however, that 7 
decreases with P for the SPC/E model, while 7 increases 
with P [Fig. ^ . This disagreement underscores the need 
to improve the dynamic properties of water models, most 
of which already provide an adequate account of static 
properties [Esf. 



VII. INTERMEDIATE SCATTERING FUNCTION 



We plo t the intermediate scattering function F{qQ,t) 



in Fig. |l^(a) for all T along the p — 1.00 g/cm'^ isochore, 
where go = 18.55 nm~^, the approximate value of the 
first peak of S{q) where the relaxation of F{q,t) is slow- 
est. We define the relaxation time r hyF{t — t) — e~^. 
We show r along isotherms in Fig. ^b), from which 
it is obvious that r has very similar behavior to D^^. 
Indeed, MCT predicts that the product Dt is constant 
along isochores, which we test in Fig. ||(c). We find 
that Dt increases slightly on cooling, but remains rela- 
tively constant along each isochore. The weak residual 
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FIG. 13. The intermediate scattering function F(q, t) (a) 
for 190 < r < 300 K along the p = 1.0 g/cm^ isochore and 
(b) for many q values at T = 210 K and p = 1.00 g/cm^. The 
solid line shows the fit to Eq. (H) for t > 2 ps. 



T-dependence in Dt should be subjected to a deeper 
scrutiny to find out if it is related to a g— vector depen- 
dent correction to scaling (since Z? is a g = quantity) or 
to the progressive breakdown of the validity of the ideal 
MCT on approaching Tc- 

The study of the time dependence of F{qo,t) allows 
us to test the predicted relation between the exponents 
b and 7 (see Eqs. (||) and (^). Since the value of b is 
completely determined by the value of 7 j26), calcula- 
tion of these exponents for SPC/E determines if MCT 
is consistent with our results. The range of validity of 
the van Schweidler power law [Eq. (^] is strongly q- 
dependent Q , making unambiguous calculation of b dif- 
ficult. 

Fortunately, according to MCT p3|, at large g- vectors, 
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FIG. 14. Fit of the stretched exponential of Eq. ^ for 
t > 2 ps at T = 210 K to both Fs^if(q,t) (o) and F{q,t) (□) 
to obtain /3. The horizontal line indicates the value predicted 
by MCT for 6 using 7 values extrapolated from Fig. [12. For 
P > 80 MPa, the relaxation of F{q, t) for g > 60 nm^^omes 
almost entirely from the first decay region, so the /? values 
obtained are not reUable in this range. 



the liquid changes significantly under increased pressure. 
To highlight the effect of structural changes on dynamic 
properties, we consider an approximately isochronic path 
- along which D remains nearly constant - such that the 
changes in dynamic properties we observe on increasing 
P are confined to their g-vector dependence. We select 
5 state points with D = (0.30 ± 0.09) x lO^G cm^/s: 
(i) T = 220 K, p = 1.00 g/cm3, (n) T = 210 K, 
p = 1.05 g/cm^, (iii) T = 210 K, p = 1.10 g/cm^, 
(i) T = 210 K, p = 1.20 g/cm3, (iv) T = 220 K, 
p = 1.30 g/cm^, and (v) T = 240 K, p = 1.40 g/cm^. 



We show in Fig. IX the g-dependence of the a-relaxation 
time T{q) extracted from the fit of F{q, t) to the stretched 
exponential of Eq. (^). For all state points, the q- 
dependence of r follows the g-dependence of S{q), as 
commonly observed in supercooled liquids and in solu- 
tions of the full (/-vector dependent mode coupling equa- 
tions. We also note that r(g) is well described by the 
relation 



T{q) cx S{q)/q^ 



(12) 



(the de Geimes narrowing relation) , as shown in the same 
figure. The MCT prediction for the (/-dependence of r is 
often very close to the relation (|l2|). 



IX. DISCUSSION 



the stretching exponent P{q), which characterizes the the 
long-time behavior of F{q,t) (see Eq. (||)), is controlled 
by the same exponent b at large q. Fits of F{qo, t) accord- 
ing to Eq. (^ are shown for many q values at T = 210 K 
and p — 1.00 g/cm^. The same fit quality is observed for 
all other low T state points. The (/-dependence of B(q) 
for p < 1.30 g/cm^ and T = 210 K is shown in Fig. |l| g|] 
for F{q,t). In addition, we show the expected value of 
b according to MCT, using the values of 7 extrapolated 
from Fig. 12, The large-g limit of (3 appears to approach 



the value predicted by MCT. Hence we conclude that the 
dynamic behavior of the SPC/E potential in the pres- 
sure range we study is consistent with slowing down as 
described by MCT [Fig. |). We also checked that the 
values of b calculated from Eq. (||) are consistent with 
the von Schweidler power law Eq. (Q), but that correc- 
tions to scaling in P'^ are relevant at several q vectors, as 
discussed in Ref. [pi. 



VIII. RELATIONSHIP OF STRUCTURE TO 
DYNAMICS 

The results shown in Fig. |l^, and the observed power- 
law dependence of diffusivity, suggest that MCT is able to 
predict the dynamical behavior of SPC/E water in a wide 
range of P and T. As discussed above, the structure of 



We have presented extensive simulations that provide 
evidence for interpreting the dynamics of the SPC/E po- 
tential in the framework of MCT. Our calculations also 
provide a necessary test of the relation predicted between 
the diffusivity exponent 7 and the von Schweidler expo- 
nent b for a wide range of values 7 and b. Our results 
support interpretation of the experimental locus Tc{P) 
as the locus of MCT transitions. 

We found that on increasing pressure, the values of 
the exponents become closer to those for hard-sphere 
(7 — 2.58 and 6 — 0.545) and Lennard- Jones (7 = 2.37 
and b — 0.617) systems thereby confirming that the 
hydrogen-bond network is destroyed under pressure and 
that the water dynamics become closer to that of normal 
liquids, where core repulsion dominates. A significant 
result of our analysis is the demonstration that MCT is 
able to rationalize the dynamic behavior of the SPC/E 
model of water at all pressures. In doing so, MCT en- 
compasses both the behavior at low pressures, where the 
mobility is essentially controlled by the presence of strong 
energetic cages of hydrogen bonds, and at high pressures, 
where the dynamics are dominated by excluded volume 
effects. We also showed how these structural changes are 
reflected in the (/-dependence of F{q,t). 

Our results underscore the need to improve the dy- 
namic properties of potentials for realistic simulations of 
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FIG. 15. Fitting parameter T(g) of Eq. (|^ in relation to S{q) along approximate isochrones. The heavy line indicates the 
prediction r ~ S{q)/q^ . 



water and other other materials. Of the many potentials 
available for studying water, only the SPC/E potential is 
known to display the power-law dependence of dynamic 
properties, but even SPC/E fails to reproduce the power 
law quantitatively. A recent study of the ST2 poten- 
tial found that the T-dependence of D is consistent 
with an Arrhenius T-dependence for T ^ 300 K, cross- 
ing over to a another region of Arrhenius behavior for 
T < 275 K j22|, in contrast to the non- Arrhenius behav- 
ior observed in real water and to our interpretation based 
on MCT for T >Tc. The presence of a low-T Arrhenius 
regime in the ST2 potential might be due to activated 
processes, that are expected to dominate the dynamics 
of fragile hquids below Tj,, as we observed for the SPC/E 
potential. Hence the ST2 potential may provide an ex- 
cellent opportunity to study these activated processes on 
a smaller time scale than is typically observed for most 
fragile liquids. 

Finally, we stress that a full comparison between the- 
ory and simulation data requires a complete solution of 
the recently proposed molecular-MCT (the extension of 
MCT to molecules of arbitrary shape) |Q. A detailed 
solution of the complicated molecular-MCT equations in 
such large region of T and P values would requires com- 
putational effort beyond the present possibilities, but a 
detailed comparison between molecular-MCT and MD 
data for one selected isobar is underway 
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APPENDIX A: FINITE-SIZE EFFECTS 

Some recent work p^ , p9| indicates that significant 
finite-size effects can affect results at temperatures close 
to the MCT Tc- Most of our systems are farther than 5% 
from Tc (i.e., T/T^-l < 0.05) so that the relatively small 
size of our system should not affect our results. We simu- 
lated two independent systems of 1728 molecules at both 
T = 200 K and 190 K a.nd p = 1.00 g/cm^ to check ff sig- 
nificant finite size effects appear at low T. The results are 
shown in Table. [II. We observe no significant deviations 
from the system of 216 molecules at T = 200 K. Hence 
we believe that no strong finite-size effects are present 
in the 216 molecule system for T > 200 K. However at 
T = 190 K, the potential energy of the larger system ap- 
pears to be significantly smaller than that in the smaller 
system. We lack adequate computer resources to make a 
reliable estimate the diffusivity in the larger system, but 
simulations are continuing in order to check the possible 
finite-size effects at this temperature. 



[1] 
[2] 

[3] 

[4] 

[5] 
[6] 



W. Gotze and L. Sjogren, Rep. Prog. Phys. 55, 241 
(1992). 

See, e.g., the special issue of Transport Theory Stat. 
Phys. 24, (1995), devoted to experimental and computa- 
tional tests of MCT. 

J.L. Barrat, W. Gotze, and A. Latz, J. Phys. Condensed 
Matter Ml, 7163 (1989); P.N. Pusey and W. van Megen, 
Physica Scripta T45 261 (1992). 

W. Kob and H.C. Andersen, Phys. Rev. Lett. 73 1376, 
1994; Phys. Rev. E 51, 4626 (1995); 52, 4134 (1995) 



W. Kob, J. Horbach, and K. Binder, :ond-mat/9812095 . 
C. Bennemann, W. Paul, K. Binder, and B. Diinweg, 
Phys. Rev. E 57, 843 ( 1998); C. Benn emann, J. 



Baschnagel and W. Paul, cond-mat/9809335, submit- 
ted, to Eur. Phys. J. B; C. Bennemann, W. Paul, J. 



[7] 

[8] 



Baschnagel and K. Binder, 2ond-mat/981002C , submit- 
ted, to J. Phys. C. 

H. J. C. Berendsen, J. R. Grigera, and T. P. Stroatsma, 
J. Phys. Chem. 91, 6269 (1987). 

P. Gallo, F. Sciortino, P. Tartaglia, and S.-H. Chen, Phys. 
Rev. Lett. 76, 2730 (1996); F. Sciortino, P. Gallo, P. 



10 



[9] 

[10] 

[11] 
[12] 
[13] 
[14] 
[15] 
[16] 



[17] 
[18] 



[19] 
[20] 

[21] 

[22] 

[23] 
[24] 



[25] 

[26] 

[27] 
[28] 



Tartaglia, S.-H. Chen, Phys. Rev. E 54, 6331 (1996); S.- 
H. Chen, P. GaUo, F. Sciortino, and P. Tartagha, Ibid 

56, 4231 (1997); F. Sciortino, L. Fabbian, S.-H. Chen, 
and P. Tartagha, Ibid, 5397 (1997). 

F.X. Prielmeier, E.W. Lang, R.J. Speedy, and H.-D. 
Liidemann, Phys. Rev. Lett. 59, 1128 (1987); Ber. Bun- 
senges. Phys. Chem. 92, 1111 (1988). 
A. P. Sokolov, J. Hurst, and D. Quitmann, Phys. Rev. B 
51, 12865 (1995). 

H. Weingartner, R. Haselmeier, and M. Holz, J. Phys. 
Chem. 100, 1303 (1996). 

C. A. AngeU, in Water: A Comprehensive Treatise, 
edited by F. Franks (Plenum, New York, 1981). 
P. G. Debenedetti, Metastable Liquids (Princeton Univ. 
Press, Princeton, 1996). 

H. Kanno and C.A. AngeU, J. Chem. Phys. 70, 4008 
(1979). 

R. J. Speedy and C. A. AngeU, J. Chem. Phys. 65, 851 
(1976); R. J. Speedy, J. Chem. Phys. 86, 892 (1982). 
P. H. Poole, F. Sciortino, U. Essmann, and H. E. Stanley, 
Nature 360, 324 (1992); Phys. Rev. E 48, 3799 (1993); 
F. Sciortino, P.H. Poole, U. Essmann, and H.E. Stanley, 
Ibid 55, 727 (1997); S. Harrington, R. Zhang, P.H. Poole, 
F. Sciortino, and H.E. Stanley, Phys. Rev. Lett. 78, 2409 
(1997). 

O. Mishima, J. Chem. Phys. 100, 5910 (1994). 

C. J. Roberts, A. Z. Panagiotopoulos, and P. G. 
Debenedetti, Phys. Rev. Lett. 97, 4386 (1996); C. J. 
Roberts and P. G. Debenedetti, J. Chem. Phys. 105, 658 

(1996) . 

M.-C. Bellissent-Funel, Europhys. Lett. 42, 161 (1998); 
O. Mishima and H. E. Stanley, Nature 392, 192 (1998); 
396, 329 (1998). 

S. Sastry, F. Sciortino, P. G. Debenedetti, and H. E. 
Stanley, Phys. Rev. E 53, 6144 (1996); L. P. N. Rebelo, 
P. G. Debenedetti, and S. Sastry, J. Chem. Phys. 109, 
626 (1998); H. E. Stanley and J. Teixeira, J. Chem. Phys. 
73, 3404 (1980). 

F.W. Starr, S. Harrington, F. Sciortino, and H.E. Stan- 
ley, Phys. Rev. Lett. 82, 3629 (1999). 

D. Paschek and A. Geiger, J. Phys. Chem. B 103, 4139 
(1999). 

H. Tanaka, J. Chem. Phys. 105, 5099 (1996). 

R. SchiUing and T. Scheidsteger, Phys. Rev. E 56, 2932 

(1997) ; T. Franosch, M. Fuchs, W. Gotze, M. R. Mayr, 
and A. P. Singh, Ibid, 5659 (1997); L. Fabbian, F. 
Sciortino, F. Thiery, and P. Tartaglia, Phys. Rev. E 

57, 1485 (1998); L. Fabbian, A Latz, R. Schilling , F. 
Sciortino, P. Tartaglia, C. Theis, |cond-mat/9812363[ 



In this work we study the motion of oxygen atoms, and 
not the center of mass motion. However, for the SPC/E 
model, the center of mass is offset from the oxygen along 
the H-O-H bisector by only 0.06 A so the difference is 
negligible. 

W. Gotze, in Liquids, Freezing, and Glass Transition, 
Proc. les Houches, edited by J. P. Hansen, D. Levesque, 
and J. Zinn- Justin (North-Holland, Amsterdam, 1991). 
F. Sciortino and P. Tartaglia, Phys. Rev. Lett. 78, 2385 
(1998). 

W. Kob, C. Donati, S.J. Plimpton, and S.C. Glotzer, 



[29 
[30 
[31 

[32; 

[33 

[34' 
[35^ 
[36 
[37 

[38' 
[39 

[4o; 

[41 

[42' 
[43^ 

[44] 



Phys. Rev. Lett. 79, 2827 (1997); C. Donati, J.F. Dou- 
glas, W. Kob, S.J. Plimpton, P.H. Poole and S.C 
Glotzer, Phys. Rev. Lett. 80, 2338 (1998). 
F. Sciortino , S. Sastry, and P. Tartaglia, cond- 



mat/9805040 



T.B. Schr0der, S. Sastry, J.C. Dyre, and S.C. Glotzer, 



cond-mat/9901271 



[45] 
[46] 



[47] 



M. Goldstein, J. Chem. Phys. 51, 3728 (1969). 

C.A. AngeU, Science 267, 1924 (1995); F. H. Stillinger, 

Ibid, 1935 (1995). 

J. P. Hansen and I. R. McDonald, Theory of Simple Liq- 
uids (Academic Press, London, 1986). 
M. Fuchs, J. Non-Cryst. Solids 172, 241 (1994). 
H. J. C. Berendsen et al, J. Chem. Phys. 81, 3684 (1984). 
O. Steinhauser, Mol. Phys. 45, 335 (1982). 
J. -P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen, J. 
Comput. Phys. 23, 327 (1977). 

L. Baez and P. Clancy, J. Chem. Phys. 101, 8937 (1994). 
S. Harrington, P.H. Poole, F. Sciortino, and H.E. Stanley, 
J. Chem. Phys. 107, 7443 (1997). 

K. Bagchi, S. Balasubramanian, M. Klein, J. Chem. 
Phys. 107, 8561 (1997). 

F.W. Starr, M.-C. Bellissent-Funel, and H.E. Stanley, 
Phys. Rev. E 60, 1084 (1999). 
F. Sciortino, unpublished. 

F.W. Starr, J.K. Nielsen, and H.E. Stanley, Phys. Rev. 
Lett. 82, 2294 (1999). 

We note that the data for D may also be fit by the Vogel- 
Fulcher-Tammann form with an "ideal" glass-transition 
temperature Tq. A detailed study of the ideal glass tran- 
sition in terms of the "configurational" entropy of the 
liquid is currently underway (A. Scala, F.W. Starr, E. 
La Nave, F. Sciortino, and H.E. Stanley, in preparation). 
M. Goldstein, J. Chem. Phys. 51, 3728 (1969). 
We estimate Tg by extrapolating Eq. ( [ll[ ) to the temper- 
ature where D reaches 10~^* cm^/s, a typical value of D 
at the glass transition. 

C.A. AngeU, J. Phys. Chem. 97, 6339 (1993); K. Ito, C.T. 
Moynihan, and C.A. AngeU, Nature 398, 492 (1999); 
F. W;. Starr, C.A. Ang eU, R.J. Speedy, and H.E. Stan- 



ley, ;ond-mat/9903451 



[49] 



[50] 

[51] 
[52] 



[53] 
[54] 



Si02, another network-forming fluid, confirms the sensi- 
tivity of the dynamics on the model potential (M. Hem- 
mati and C.A. AngeU, in Physics meet Geology, edited by 
H. Aoki and R. Hemley (Cambridge Univ. Press, Cam- 
bridge, 1998)). 

We do not show the behavior of /3(g) for p — 1.40 g/cw? 
and T — 210 K since the MCT predictions are not ex- 
pected be valid for this state point, as discussed in Sec. ^ 
T. Franosch, M. Fuchs, W. Gotze, M.R. Mayr, and A. P. 
Singh, Phys. Rev. E 55, 7153 (1997). 
P.G. de Gennes, Physica 25, 825 (1959). 
For hard spheres, see J.L. Barrat, W. Gotze, and A. 
Latz, J. Phys. Condensed Matter Ml, 7163 (1989); For 
Lennard- Jones, see U. Bengtzelius, Phys. Rev. A 34, 
5059 (1986). 

F.H. Stillinger and A. Rahman, J. Chem. Phys. 60, 1545 
(1974). 

L. Fabbian, A. Latz, R. SchiUing, F. Sciortino, P. 
Tartaglia, C. Theis, Phys. Rev. E., in press. 



11 



TABLE I. Summary of the state points simulated with 216 molecules interacting via the SPC/E potential. For all state 
points, the uncertainty in the potential energy U is less than 0.05 kj/mol. The uncertainty in the diffusion constant D is 
approximately ±4, in the last digit shown. State points we equilibrated for a time teq, followed by "data collection" runs of 
duration (jata- 
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o 
o 




I.oU 


cn 1 /I 
— 0U.i4 


1 1 OT _1_ 1 /I 

llzr it 14 


Q C/1 

O.04 


c 



o 
o 




1 /in 


/i n 07 

— 4y.y ( 


1 n'yn -L i /i 
ly / y ± 14 


1 Qn 


c 



o 
o 


300 


0.95 


—46.80 


— 109 ± 12 


19.9 


0.5 






1.00 


—47.20 


—13 ± 13 


20.0 


0.5 


j 




1.05 


—47.49 


112 ± 14 


18.3 


0.5 






1.10 


-47.65 


264 ± 14 


18.2 


0.5 






1.20 


-47.95 


678 ± 16 


15.3 


0.5 






1.30 


-48.06 


1293 ± 18 


11.2 


0.5 






i.lO 


-i7.cS8 


2222 ± 11) 


1.95 


0.5 




350 


0.90 


-43.21 


-105 ± 16 


61.1 


0.5 


40 ps 




1.00 


-44.35 


62 ± 18 


49.7 


0.5 


40 ps 




1.10 


-45.15 


358 ± 20 


38.1 


0.5 


40 ps 




1.20 


-45.56 


828 ± 22 


27.0 


0.5 


40 ps 




1.30 


-45.76 


1504 ± 25 


18.0 


0.5 


40 ps 




1.40 


-45.50 


2522 ± 26 


13.9 


0.5 


40 ps 
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TABLE II. Fitting parameters to the power law predicted by MOT, for T < 300 K. Note that the state points 
{p = 1.00 g/cm ^ T = 200 K and T = 190 K) and {p = 1.40 g/cm ^, T = 210 K) are not included in the fit because 

they arc very close to Tc, and so do not conform to the power law. 



p (g/cm^) or P (MPa) 


Do (10-® cmVs) 


TciK) 


7 


0.95 


159 


201.4 


2.84 


1.00 


114 


193.6 


2.80 


1.05 


78.5 


188.3 


2.67 


1.10 


63.9 


188.6 


2.31 


1.20 


54.0 


189.9 


2.24 


1.30 


57.3 


192.9 


2.70 


1.40 


47.8 


210.2 


2.59 


-80 


133 


197.9 


2.79 





102 


193.9 


2.62 


100 


77.2 


187.8 


2.61 


200 


65.1 


188.8 


2.50 


300 


60.8 


188.8 


2.42 


400 


59.1 


190.7 


2.25 



TABLE III. Summary of the state points simulated with 1728 molecules to check for finite-size effects. For both state points, 
the uncertainty in the potential energy U is less the 0.04 kj/mol. 



T 


P (g/cm^) 


U (kJ/mol) 


P (MPa) 


D (10-® cmVs) 


teq (ns) 


<data (ns) 


190 


1.00 


-55.12 


6±8 




60 


40 


200 


1.00 


-54.41 


-4 ±7 


0.0020 ± 0.0007 


30 


35 
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